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We propose to use quantum tomography to characterize the state of a perturbed Bose-Einstein 
condensate. We assume knowledge of the number of particles in the zero-wave number mode and 
of density distributions in space at different times, and we treat the condensate in the Bogoliubov 
approximation. For states that can be treated with the Gross-Pitaevskii equation, we find that 
the reconstructed density operator gives excellent predictions of the second moments of the atomic 
creation- and annihilation operators, including the one-body density matrix. Additional inclusion of 
the momentum distribution at one point of time enables somewhat reliable predictions to be made 
for the second moments for mixed states, making it possible to distinguish between coherent and 
thermal perturbations of the condensate. Finally, we find that with observation of the zero-wave 
number mode's anomalous second moment, (fio&o + "o^J), the reconstructed density operator gives 
reliable predictions of the second moments of locally amplitude squeezed states. 
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I. INTRODUCTION 

Bose-Einstein condensates offer countless demonstra- 
tions of macroscopic quantum effects pj, Q. Recently, 
there has been growing interest in the use of condensates 
as quantum optical components and as atom laser beam 
sources. The effective use of very many of these applica- 
tions requires knowledge of certain aspects of the conden- 
sate's quantum state, and several articles have dealt with 
theoretical calculations of the states obtained by various 
perturbations of a condensate. 

We present here a method by which experimental data 
can be used to estimate characteristics of the system's 
quantum state. For this purpose we generalize the ap- 
plication of the Jaynes principle of maximum entropy 
(MAXENT) [3j for quantum state tomography, aspro- 
posed in Q and applied to single atoms in [j| and p, to 
the dynamics of a condensate described in the Bogoliubov 
approximation. We assume that the fraction of particles 
in the zero-wave number mode and the position distri- 
butions at different times after preparation of the sys- 
tem have been measured. It is important to note that 
we are not claiming that the MAXENT density opera- 
tor is an approximation to the true many-body operator 
of the quantum state; indeed this would surely require 
more than just density measurements (one-body opera- 
tors). Instead, we seek to correctly predict the second 
moments of the atomic creation- and annihilation oper- 
ators 4>(x) and tp^x). Because these second moments of 
the ladder operators will be discussed extensively below, 
we shall call second moments of the form [rj}< (x)ijr (a/)) 
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and (ip{x)ip(x')) anomalous second moments, while sec- 
ond moments of the form (x)xjj(x')) will be referred 
to as normal second moments. 

As an example of applying the method, three types 
of quantum mechanically very different states will be 
studied. First, we study states characterized by a sin- 
gle Gross-Pitaevskii wave function with small deviations 
from the constant value attained by the ground state 
homogeneous condensate. Then secondly, we apply our 
method to an incoherently exited system, modeled by the 
addition of a localized thermal component to the homo- 
geneous condensate. Finally, as a third case, we study the 
method's fidelity for a system with a localized amplitude 
squeezed atomic field. 

Our paper is organized as follows. In section[n]we give 
a brief introduction to the central ideas behind tomogra- 
phy and a description of the procedure for reconstruct- 
ing the density operator of a quantum system by use of 
Jaynes principle of maximum entropy. In section ITTT1 we 
specify our assumptions about the condensate, and we 
introduce the Bogoliubov approximation and the general 
para-unitary transformation procedure. In section llVl we 
will combine the MAXENT and Bogoliubov theories. In 
section we apply the machinery to the three different 
trial states. In section IVTl we discuss the applicability of 
the method and we finally conclude the paper in section 

eh 



II. QUANTUM STATE TOMOGRAPHY 

A. Tomography 

Can you find all the peaks of a mountain landscape or 
draw a map of a city just by seeing its skylines from all 
compass directions? Of course not: obstacles may block 
your line of vision, no matter from which angle you look, 
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making it impossible to determine the layout. The cen- 
tral idea behind tomography is the, somewhat surprising, 
fact that if you had performed absorption measurements 
from all angles, instead of looking at shadows, you could 
indeed have found every detail of the area. 

Finding its first application in medicine, tomography 
techniques have been used extensively in this field since 
the 1970's. Here, tomography is used to reconstruct 3-D 
pictures from data recorded using X-rays or NMR. 

With regard to quantum systems, tomography can be 
used to reconstruct the density matrix of the system. For 
quantum harmonic oscillator systems this reconstruction 
can be done exactly by measuring the quadrature dis- 
tributions for all quadrature angles; for example mea- 
suring the density distribution in space for all points of 
time during one period of oscillation. The techniques 
for this include the inverse Radon transformation Q, 
and the more reliable technique of pattern functions 
The analogy to the absorption measurements mentioned 
earlier is apparent if one considers the quadrature distri- 
butions as marginal distributions of the system's Wigner 
function. 

Buzek and Drobny have shown that a good approxima- 
tion to the density matrix for the one-particle harmonic 
oscillator can be achieved using the MAXENT princi- 
ple with the observations of just a few quadrature dis- 
tributions and the mean oscillator excitation number £| . 
Skovsen et al. recently used MAXENT tomography to 
reconstruct the one-body density matrix of free particles 
from experimental data |6(]. 

In this paper we will extend the scope of MAXENT 
quantum tomography by using it to approximately re- 
construct the second moments of the ladder operators, in- 
cluding the one-body density matrix, of a many-particle 
system: a condensate treated in the Bogoliubov approx- 
imation. 



B. MAXENT principle for reconstruction of 
density operators 

The method we use to estimate the quantum state of 
our system at hand is by finding the MAXENT den- 
sity operator as suggested by Jaynes 3]. The goal of 
the method is to find a density operator that exactly re- 
produces the knowledge one has about the system while 
making as few assumptions as possible about degrees of 
freedom one has no knowledge about. A proposal for this 
is the maximum entropy density operator pme chosen for 
having maximal von Neumann entropy S = —Tr(phip). 
We make this more rigorous in the following way |2*cj : 

Let a set of observables {G v }, (y £ 1, 2, . . . , n) be asso- 
ciated with a quantum system prepared in an unknown 
state p. We imagine an ensemble of these quantum sys- 
tems wherein there are no correlations between the sub- 
systems of the ensemble, and we assume to have mea- 
sured the ensemble expectation values {G v } of the ob- 
servables {G„}- 



Unless the set of observables {G u } constitutes the quo- 
rum (a complete set of observables) there will in general 
be many density operators P{q^ that satisfy the normal- 
ization condition Tr(/3^j) = 1 and predict the measured 
values: 

Tr(p 0} G v ) = G v , ^€{l,2,...,n}. (1) 

The MAXENT principle says that the most unbiased 
guess for the density operator approximating p and ful- 
filling these conditions is the one with maximal von Neu- 
mann entropy: 

S(pme) = max[S{p {d} )-yp {d} ], (2) 
S(p) = -Tr(phxp). (3) 

As shown by Jaynes this implies a density operator of 
the form: 



Pme 



j— exp ( - \ V G V ) 

J {G} \ u=l J 



(4) 



with the generalized partition function: 



Z 



{G} 



Tr 



exp 



v=l 



\ V G V 



(5) 



and with the {Aj,} being a set of Lagrange multipliers 
chosen so that pus in Eq. (QJ fulfills the conditions 
Eqs. CEJ. Using these A„'s we have the explicit form of 
the MAXENT density operator Eq. gj). 

A technical difficulty arises because, in all but the most 
trivial circumstances, it is not feasible to analytically in- 
vert the conditions Eqs. given Eq. (0J to obtain the 
set {A„}. For implementation it is more convenient, and 
we shall indeed use this procedure, to minimize the norm 
of the differences between measured values and values 
predicted by pme with respect to {A^}: 



/ Tr(p ME ({\ v })G l )-G l 
\Tr(p ME ({K})G n ) - G n 



(0) 



Here, we have given the same weight to all components 
of the vector of errors. One may wish to modify this if, 
for instance, some observations are made with a much 
better precision than others. 



III. DESCRIPTION OF THE CONDENSATE 

A. Bogoliubov Hamiltonian 

Let us consider a 3-dimensional cold gas of atoms hav- 
ing mass m, occupying a volume V, and dominated by 



3 



s-wave collisions where U(r,r') — UoS(r — r') is the inter- 
particle potential. Such a gas is well represented by the 
Hamiltonian: 



H 



full 



e k a\a k + - ^ a^+X'-^*'®* (7) 



k,k' ,q 



where 



is the free-particle energy and g 



Uo 



the coupling constant. The operators a k and a k are the 
usual boson ladder operators annihilating and creating a 
particle with wave vector k, and having the commutation 
relations: 



[a k ,a\,] = Sk,k' 
[&k,ak r ] = 0. 



(8) 
(9) 



A Hamiltonian similar to Eq. Q applies in one and two 
spatial dimensions with suitable redefinitions of g, de- 
pending on the confinement of the remaining coordinates 
|lO| , |lf| . For simplicity, only one spatial dimension will 
be used in this paper. The wave vectors k are then merely 
scalars, and will in the following be referred to as wave 
numbers. 

Due to the complicated dynamics of the Hamiltonian 
Eq. (JZJ), it is often necessary to use an approximate 
Hamiltonian to make analytical calculations possible. 
We will employ the so-called Bogoliubov approximation, 
which is valid when the vast majority of the particles are 
in the zero wave number mode [l2j: 



(aja ) > ^(a£a fc ) 



Having this condition fulfilled, it is a good approximation 
to approximate the Hamiltonian Eq. {JJ by |21| : 

H = - ^2 ( (e k + gNtot) (a\a k + a k a f k ) + 

fe^O ^ 

+ gN tot (a k a- k + al fe al)^ (11) 

and Ntot is the total number of particles in the gas. This 
operator is bilinear in the ladder operators, and can be 
diagonalized by a change of basis. We will do this by 
using elementary linear algebra. 

B. Discretization of position and momentum 

We will now specify some useful technical details of 
our treatment. We consider an odd number (N > 1) of 
evenly spaced points in 1-D coordinate space, symmet- 
rically distributed around (and including) the origin, on 
an interval of length L = NAx. 



-M,-M+1,...,M (12) 



where M = ^j^. We assume that the spatial atomic 
density is measured on this grid, i.e., our observations 



will be of the form (tf(x n , t)ijj(x ni t)). Due to the form 
of the Hamiltonian Eq. , it is convenient to introduce 
the discrete wave numbers corresponding to the spatial 
discretization above. 



2vr 

k q — qAk, Ak — — , 
L 



(13) 



where q can assume the same values as n above. 

The usual discrete Fourier transform from coordi- 
nate space (i/j(x n ) and ^(xn)) to wave number space 
(a(kg) and a'(k q )) reads: 



a(k q ) = — L= ^2 ip(x n ) exp (ik q x n ) (14) 

^ n=-M 
1 M 

a)(k q ) = — = ^ ^(x n )e^q>{-ik q x n ). (15) 

^ n=-M 



All the ladder operators in coordinate space can be 
arranged in a 2N x 1 column vector: 



4> 



4>(x M ) 

^{X-m) 



(16) 



(10) 

and correspondingly: 



In this paper, we shall denote vectors of this type by bold 
letters. Similarly, we can construct the column vector 
a from the operators in wave number space a{k q ) and 



aHk q y. 



( a{k- M ) \ 



a(k M 
at(fc 



(18) 



V a\k M ) J 

Using these vectors, we can write the Fourier transfor- 
mation Eqs. H14 |l - 115(l as a matrix multiplication by a 
2N x 27V matrix: 

a = ■ ip, whereby a' = xjy ■ (19) 

Here, and in the following, we shall use script letters to 
denote matrices of dimension 2N x 2N. 
The matrix si has the following form: 



a/ = 



F 
F* 



(20) 
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and F is the usual N x N discrete Fourier transformation 
matrix, whose g'th row is: 



v N 



M ),---,exp(ik q x M )] (21) 



from which Eqs. (|14fl - i|15[) are easily recovered. 

It is straightforward to see that the operators a(k q ) and 
ip(x n ) obey similar commutation relations Eqs. JSJ - ©- 
In this paper we shall exclusively deal with transforma- 
tions to new sets of operators that conserve this property. 



C. Diagonalization of the Bogoliubov Hamiltonian 

Using the notation developed in the preceding section, 
we can write the Bogoliubov Hamiltonian Eq. in a 
compact form: 



(22) 



As shown in e.g. [13|], this Hamiltonian can be diago- 
nalized to the matrix $ by transforming to new boson 
operators: 

b = 38a (23) 



By diagonalizing the Hamiltonian we have found a ba- 
sis of non-interacting modes with a simple time evolution: 



b q (t) = b q (0)exp(-ie(q)t/H). 

By defining the para- identity matrix ^ , which will play 
an important role in the next section: 



-1) 



(27) 



diag(l, ...,1,-1,. 

N N 

we can write the vector of operators b at time t as: 

b{t) = exp (-iJn/K)b(0) = W(t)b(0). (28) 
The changes of basis can be summarized like this: 
■0(f) = ^" 1 a(*) 



^- 1 ^~ 1 'W(t)b(0), 



(29) 



whereby we can find the time evolution of the tp(x n ,t) 
and f (i n , t): 



(30) 



- „t 



b\&) 
b^b, 



M'SS^^b 



(24) 



(25) 



The matrix J(f only has elements different from zero in 
the two diagonals Jf?(k q , k q ) and J4?(k q , k- q ), giving the 
matrix an " X-structure" . In addition, all elements are 
real. Therefore the matrix SS~ X that diagonalizes 3^ also 
only has elements different from zero with these indices, 
and can be chosen real. The upper-left to lower-right 
diagonal elements are denoted u(k q ), while the upper- 
right to lower-left diagonal elements are denoted v(k q ). 
The elements can be found in e.g. |l4j: 



1 



U(k q ) = < / g 2 /4+7/2 



<q)/E 



V{k q 







for q = 
for q 7^ 



for q = 



£74+7/2 



i for q 0. 

The diagonalized form § has elements t{k q ): 



where E = , 7 



e(k q ) = Eoy/ryqZ + cf/4: (26) 
gNtot , and gN to t has the same 



meaning as in Eq. (|llfl . 



-Bo 



D. General para-unitary diagonalization 

All the vectors of operators Eq. JT5J, Eq. i|T%)l and 
Eq. lt2*3l) are examples of a vector of general boson ladder 
operators: 



a 



olm 

t 

M 



(31) 



fulfilling the commutation relations: 



[a q ,a q >] = 



(32) 
(33) 



Both the Fourier transformation Eq. I|19() and the 
Bogoliubov transformation Eq. H23J) ar< 3 examples of 
changes of basis, known as para-unitary transformations, 
which conserve the Boson commutation relations between 
the operators. All the transformations to new sets of 
operators in this paper conserve this property. In our 
present treatment, there are two main reasons for mak- 
ing these transformations. Firstly, the measurements we 
have performed on the system are in the coordinate space 
and wave number space at different times, but the oper- 
ators corresponding to these measurements have compli- 
cated time evolutions. Working instead in the 6-basis 
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Eq. (|28|l the time evolution is simple, and the exponent 
in Eq. (0J can be easily calculated. Secondly, to per- 
form the traces in Eq. it is convenient to shift to 
a basis where the density operator is diagonal, and this 
can be done by a para-unitary transformation. We recall 
that there are 2N ladder operators, regardless of basis 
(f.x. the set {tp(x n )} and {i/j'(x n )}), which will make 
the transformation matrices 2N x 2N. 

The commutation relations demand that in a transfor- 
mation from one set of boson operators to another: 

(3 = ,9 a. (34) 

the 2N x 27V matrix 3" must be para-unitary, i.e. satisfy 
the condition |22l ]: 

J = (35) 

with the diagonal matrix defined in Eq. (|27|l [23| . 
From the definition it is seen that the product of two 
para-unitary matrices is again para-unitary. 

The matrices we will diagonalize in this paper will all 
be hermitian, positive definite and have the structure 

& = ( q* p* ^ • (36) 

Let X be a hermitian operator in the a. basis with co- 
efficient matrix ' . Such a hermitian matrix X can be 
diagonalized by the para-unitary matrix 3 to a new set 
of operators {/3 ra } and {/3j }: 

X = ctXa (37) 

= pi&p. (38) 

Because JT is positive definite the para-eigenvalues Jzf^i 
will also be positive by Sylvester's law of inertia. The 
matrix 2? can be chosen to have the structure: 

*=(vr) (39) 

giving Jzf the form: 
_S? = diag{££\ > \,££'ifr . ■ ■ , JSfjv.jv, JSfi,i, ■ ■ ■ ,J£n,n)- (40) 

This general procedure will become useful in the fol- 
lowing section, where we will need to evaluate expecta- 
tion values of the observed densities within a quantum 
state of the form Eq. facilitated by a para- unitary 
diagonalization of the exponent ^ v \ V G V . 

IV. USING MAXENT WITH THE 
BOGOLIUBOV APPROXIMATION 

Now we are ready to use the MAXENT formalism on 
the condensate in the Bogoliubov approximation. As the 



set of observables {G v } we use measurements of particle 
numbers at coordinate space points at different times t, 
and of the occupation of wave number modes at times if . 
The times t can, but need not, be equal to the times t' 

{G v } = {^(x n ,t)^(x n ,t)+^(x n ,t)^(x n ,t)}[j 

{a){k q , t')a{k q , + a{k q ,t')a\k q ,t')}. (41) 

The reason for using this symmetrical form of the observ- 
ables is that it leads to the favorable structure Eq. i|36|) 
of the matrix in the exponent of Eq. (0J . 

The task is now to put these observables into the 
MAXENT density operator Eq. © and perform traces 
like Eqs. 0. To do this we will use the same method as 
in section ^ to diagonalize the matrix in the exponent 
into a set of new, non-interacting quasi-particles. 

Let W n be the 2N x 2N matrix with the following 
elements: 

{1 forr = s = n + M + l 
1 for r = s = n + N + M + l (42) 
otherwise 

so that (using Eq. (|29|l for the first part): 

ft{x n ,t)$[x nt t) + 4>(x n , t)ft(x n , t) = 4>\t)W n iP(t) 



= &t(o)«rt(t) (^- 1 ) t (s/- i yw n */- 1 a- 1 &(t)b(o) 

(43) 



a 1 (k q ,t')a(k q ,t') + a{k q ,t')a} \k q ,t') = a f {t')W q a(t') 



= &t(0)* rt (t / ) W q ®- l( %{t')b{Q). (44) 

The operator ^2 V \ V G V , with G v being operators of 
the form Eq. (|43|l and Eq. (|44|l is hence expressed in 
terms of the Bogoliubov eigenmode operators, and we 
can write the MAXENT density operator Eq. J3J in a 
compact form: 

Pme = iexp(-&t(0)5»&(0)), (45) 

where 

M 

& = E E a m- 

t n=-M 

w\t) (^- 1 ) t (^/^ 1 ) t w n $t- x m- x w{t) 

M 

+E E 

t> q=—M 

W\t') (^- 1 ) f l^ q ^- 1 ^(t'). (46) 
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Since the structures of *% (t) , $4 and 38 are all like that of 
Eq. and the structure of W q is like that of Eq. (gSJ, 
so is the structure of & . We can therefore diagonalize 
it to a diagonal matrix Jz? by a para-unitary change of 
basis to new Boson operators c = ^6(0): 



6^0)^6(0) 



N 



j + c i c jJ 



(47) 
(48) 



Remembering that the are all positive, the partition 
function Eq. (J5J becomes: 



Tr exp 



n 



c 7 H~ Cj Cj- 



L Vexp^j) - exp(-_Sfj 



(49) 



and we can readily determine the corresponding set of 
expectation values: 



(c m c rl ) — 



\ Tr I c] n c n exp 



dhxZ 



N 



(Cm Cn) 



exp(2Jz? n ,„) - 1) 



= (44) = o. 



and 



(50) 



(51) 



Knowing these expectation values, it is easy to find the 
expectation values of all the second moments in any basis, 
for example: 



here that we must require the diagonal elements of jSf to 
be strictly positive for the geometrical series summed to 
find Eq. (|49|) to be convergent. Therefore, by Sylvester's 
law of inertia, 2? must also be positive definite. In ad- 
dition, for this condition to be fulfilled, the number of 
observables in {G„} must be at least as large as N. The 
reason for this is that in Eq. I|4fci[l the matrices Wn have 
dimensionality two in the 2N space, the dimensionality of 
course being conserved under the para-unitary changes of 
basis. So to span the 27V space, at least N observables of 
dimensionality two must be included. In the present pa- 
per, however, we will have many more observations than 
N, and linear dependencies between them are highly un- 
likely to reduce the dimensionality below 2N. 



TOMOGRAPHY ON THREE QUANTUM 
STATES 



A. Initial excited states 

The time evolution is given by Eq. i|3U|) . Since we 
are not attempting to reconstruct the full many-particle 
quantum state, but only the second moments of the lad- 
der operators, it is hence sufficient to specify these at, for 
example, t = 0: 



^ 1 ^ 1 ^(i)^^(V(0)i/' t (0)) 



(52) 



To examine the reliability of the MAXENT technique, 
regarding reconstruction of the correct second moments 
of the ladder operators, we study three very different 
types of condensate quantum states with nearly the same 
density distribution at t = 0: namely a gaussian pertur- 
bation on top of a flat condensate. 



{^(x-m^Hx^m)) ■■■ (-0 t (x_ M )'0(a;-M)) • 



V 



■••/ 



= ss?- 1 {a{t)a) (i)) (^" 1 ) t 
= [^-\t)3g£/y 1 (cc^) ■ 

In passing, it is noted that the lower right N x TV sub- 
matrix of the matrix {^(t)^ (t)) is the one-body density 
matrix. 

As a final point we give a comment on the require- 
ment of the positive definiteness of & . It is clearly seen 



1. Coherent state is a pure state of all the parti- 
cles, representable as a product state with all the 
particles in the same state. We choose this to be 
a constant with a gaussian perturbation in coordi- 
nate space. 



\Statel) 

4>l 



(4 

E 



4 |0), where 

gauss {Xn )^P^ (Xn ) 



(53) 
(54) 



and 



'gauss 



is a constant plus a gaussian with stan- 



dard deviation a, centered at the origin: (j} ga uss oc 

v 21 



1 + 7] exp 



(^) 2 ]} 



N, 



part 



is the total num- 



ber of particles. We can envision this state formed 
in a uniform system with a negative potential dip, 
having State 1 as its ground state. The spatial evo- 
lution of the state after abruptly turning off the 
potential at t = is what is usually handled with 
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the Gross-Pitaevskii equation: 



(a) 



-ih 



at 



h 2 d 2 

2m dx 2 



+ g N tot \^\ 2 )^ 



(55) 



2. Thermal perturbation This state is a flat con- 
densate with a thermal gaussian perturbation su- 
perposed hereupon. The gaussian perturbation is 
treated as originating from a thermal boson gas in a 
harmonic trap and we have assumed no initial cor- 
relations between the perturbation and the original 
flat condensate. Therefore the second moments for 
this state in coordinate space are simply the sum 
for the thermal gaussian and the flat condensate. 

3. Squeezed condensate We get this state by apply- 
ing an operator, similar to the squeezing operator 
known from squeezing of light, to a condensate with 
all Np ar t particles in the condensate mode: 

\State3) = exp z*(^ g ) 2 - z{^ g ) 2 (Sj) |0). (56) 

The most important difference between the two 
preceding states and this one is the magnitude of 
the anomalous second moments at t = 0. In State 1 
and 2 the anomalous second moments are all zero, 
where they in this state are comparable in magni- 
tude with the normal second moments. 



B. Results of tomography 

We present here the results of using the MAXENT 
technique on the three states discussed in section IVAl 
The main points of interest are not just whether a reliable 
reconstruction is attained, but also the amount of data 
needed for this. In all the results shown in this section 
we have used density distributions in coordinate space 
at four different times. In State 2 we have furthermore 
used the momentum distribution at t = 0, and in State 3 
we have additionally included an observation of (cioQo + 
ajaj). 

We have used 10 5 particles in the flat part of the con- 
densate at t = and set gN to t =0.1. The population in 
the zero wave number modes is 99.5% in State 1 and 3 
and 96% in State 2. For the calculations we have used a 
grid of N = 25 discrete points. 

In Fig. n we show the density distributions in coordi- 
nate space of the three different states declared in section 
IV Al We have intentionally chosen the initial distribu- 
tions to be very similar. The perturbations give rise to 
density variations of a magnitude readily detected 0] . 

We now proceed to examine the true second moments 
of the ladder operators for the three states and those 
reconstructed using the MAXENT procedure. We will 
study these using surface plots of the absolute value of 
the second moments and of the absolute value of the dif- 
ference between the true values of these moments and 




X / A X 



(b) 




t/[h/(27t*E n )] 



X / A X 



(C) 




X / A X 



Figure 1: Density distributions in coordinate space at four 
different points of time for: (a) Statel, (b) State2 and (c) 
StateS. The data shown are those used in the calculations 
together with momentum distributions at t = for State 2 
and 3, and one anomalous moment for State 3. 



those predicted by MAXENT for the state at t = 0. In 
this way, errors in both magnitude and phase of the ma- 
trix elements will be visible. 

For State 1, where we will see that the reconstruction is 
precise and only the normal second moments are non-zero 
at t = 0, we shall also display the phase-space Wigner 
function corresponding to these normal second moments. 
Apart from normalization, the normal second moments 
can be interpreted as the density matrix for any single 
atom in the system, and the phase space Wigner func- 
tion W(x n , kg) is a convenient representation of this. The 
cited works p[, all show the Wigner function of the 
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particle states, and much interest has been devoted to the 
fact that measurements of exclusively positive marginal 
distributions may be used to identify Wigner functions 
with domains of negative values ^(|, ^tJ- The connec- 
tion between the Wigner function and the normal second 
moments is given in appendix iBl 



1. State 1 

Treating first State 1, Fig.[3presents the absolute value 
of the normal second moments as found from section lV Al 
The reconstructed second moments are identical to the 
input values within roundoff errors, which amounts to 
an accumulated error of 10 on each matrix element, 
including phase. The same is true for the anomalous 
second moments, all having the true value zero at t = 0. 
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Figure 2: Absolute value of the reconstructed normal second 
moments (%}>' (x r )tf>(x s )) of Statel. The anomalous second 
moments are all zero (not shown). Both normal and anoma- 
lous moments are reconstructed to within machine precision, 
including complex phase, giving accumulated errors of magni- 
tude 10 -4 on each element. The data used for the reconstruc- 
tion are four spatial density distributions. The true second 
moments are not shown, as they are practically identical to 
the reconstructed ones. 



Figure 3: The upper graph (a) shows the Wigner function 
for the normal second moments of Statel, i.e. the one-body 
density matrix apart from normalization. In the lower graph 
(b) the elements W(n,0) are set equal to zero to more clearly 
reveal the finer structures and the negativities of the Wigner 
function. 



In Fig. [5] we display the absolute value of the recon- 
structed anomalous second moments, that all have the 
true value zero at t = 0. 



As a consequence, the Wigner function for State 1 is 
also perfectly reconstructed, and is presented in Fig. |3J 



2. State 2 

For State 2 the true and reconstructed normal second 
moments are shown in Fig. Also shown in this figure 
is the absolute value of the difference between the true 
and reconstructed normal second moments. The errors 
are many orders of magnitude larger than for State 1. 



3. State 3 

The absolute value of the true and reconstructed sec- 
ond moments are shown in Fig. [|J] together with the ab- 
solute value of the difference between these. 

In State 1 and 2, the true value of the anomalous sec- 
ond moments have all been equal to zero at t = 0. In con- 
trast, for the squeezed state, the absolute value of these 
elements are comparable with the normal second mo- 
ments. The absolute value of the true and reconstructed 
second anomalous moments are displayed in Fig. {7\ 




Figure 4: The first two graphs show the absolute value of (a) 
the true, and (b) the reconstructed normal second moments 
{$(x r )%j){x s )) for State2. The errors above the condensate 
are up to about 10% of the peak value. The bottom graph (c) 
shows the absolute value of the difference between true and 
reconstructed normal second moments. The data used in the 
reconstruction are four spatial density distributions and one 
momentum distribution. 



Figure 5: Absolute value of the reconstructed anomalous sec- 
ond moments (tp(x r )^>(x s )} for Statel at t = 0. The true 
values are all zero at t = 0. The magnitude of the errors is 
similar to the errors of the normal second moments for this 
state. 



VI. DISCUSSION 

Initially, a short comment on the choice of observables 
might be in order. If there were no constant background 
condensate the momentum distribution could in princi- 
ple be found from the spatial distributions at late times. 
Having to deal with this background, however, one has to 
include the additional observables in the MAXENT den- 
sity operator to exclude the possibility of the flat part of 
the condensate being an almost even distribution of par- 
ticles with all allowed types of wave numbers whirling 
left and right, which would otherwise be preferred by the 
MAXENT formalism. 

Indeed we find for State 1 a very good reconstruc- 
tion, the errors being of the same magnitude as com- 
puter roundoff errors, from using just measurements of 
the distribution in coordinate space at three times and 
the number of particles in the zero wave number mode 
(a^(O)a(O)). In addition to Statel in section IV Bl we 
also tried to reconstruct states of this type with more 
complicated initial distributions, but still describable by 
a simple Gross-Pitaevskii wave function. These included 
perturbations with two peaks and perturbations with one 
peak and one hole in the flat background. In all cases the 
reconstruction had the same precision as Statel, indi- 
cating that the second moments, including the one-body 
density matrix, of states describable by a simple Gross- 
Pitaevskii equation can be completely reconstructed. 

For State 2, we found that a wave number distribu- 
tion also had to be included as observable to get a some- 
what reliable reconstruction of the second moments. It 
is reasonable that it is difficult to tell apart whether the 
moving particles belong to the perturbation on top of a 
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Figure 6: The first two graphs show the absolute value of (a) 
the true, and (b) the reconstructed normal second moments 
(ffl (x r )xj)(x s )) of State?,. The lowest graph (c) is the absolute 
value of the difference between the true and reconstructed 
second moments. The data used in the reconstruction are four 
spatial density distributions, one momentum distribution and 
(oofio + a' a' ). 



Figure 7: The first two graphs show the absolute value of (a) 
the true, and (b) the reconstructed anomalous second mo- 
ments {^)(xr)i){x s )} of StateZ. The lowest graph (c) is the 
absolute value of the difference between the true and recon- 
structed second moments. The data used in the reconstruc- 
tion are four spatial density distributions, one momentum dis- 
tribution and (ao&o + a' a' ). 
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flat condensate or to some collective motion of the whole, 
for example similar to State 1. Including the momentum 
distribution remedied this problem somewhat, but errors 
still persist (see Fig.QJ. Subtracting the flat background, 
these errors are of magnitudes about 10% of the pertur- 
bation peak value. 

For State 3 with large numerical values of the anoma- 
lous second moments, we find a reasonable reconstruction 
at t = of the normal second moments, but reconstruc- 
tion of the anomalous second moments completely fails, 
unless we additionally include the observable doa,Q + ajaj 
in the MAXENT density operator, as we did in section 

EH 

A squeezed condensate may be useful in atom inter- 
ferometry as input to the dark side of an atomic mir- 
ror, whereupon a bright coherent input will be split with 
number fluctuations below the values for a binomial dis- 
tribution It is in a similar setup that (aoao + ^o^o) 
can be determined experimentally. 

It is also worthy of notice that there was almost no 
improvement in the reconstruction from including ad- 
ditional position and momentum distributions for more 
points of time, and that using also the observable (aoao + 
agaj) for State 1 and 2 did not change the precision of 
the reconstructed second moments. 

As for states significantly different from the ones 
treated here, distributions in coordinate space at more 
points of time may be required for a reliable reconstruc- 
tion. Since enlarging the number of measured distribu- 
tions in space is rarely a problem, we suggest for specific 
application to gradually include more distributions until 
the results have stabilized. The calculation time is quite 
manageable, the shown examples all taking well under an 
hour on an ordinary PC. 

VII. CONCLUSION 

We have shown how to combine the MAXENT 
principle with the Bogoliubov approximation to give 
reliable estimates of the normal and anomalous second 
moments of the ladder operators, e.g. (x)ip{x')) 
and (^(a^t/^a;')), in a perturbed condensate. For states 
describable by a simple Gross-Pitaevskii equation, we 
find near perfect reconstruction from data of density 
distributions at a few points of time and the number 
of particles in the zero wave number mode. For states 
with a larger amount of particles in non-zero wave 
number modes compared to the size of the perturbation 
in coordinate space, we find it necessary to include 
a momentum distribution to obtain a somewhat reli- 
able reconstruction. Errors in this case are less than 
10% of the perturbation peak value (peak above the 
constant condensate background). Thus, the method 
makes it possible to distinguish between a thermal and 
a coherent perturbation of the condensate. Finally, 
we find that by measurement of a single observable 
related to squeezing, apart from measurements of 



momentum and coordinate space distributions, the 
method is able to correctly reconstruct all second 
moments of the ladder operators for squeezed states. 
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Appendix A: COMPUTATIONAL DETAILS 

Several numerical problems may arise in finding the 
set {A y } that satisfies Eqs. ffl. As already touched upon 
towards the end of section lll Bl we have treated the prob- 
lem as a minimization of deviations. We have chosen to 
minimize the vector of deviations Eq. © by using the 
Levenberg-Marquardt algorithm in the Matlab package. 
The most common problem encountered when using this 
approach is that for a given trial set {A„}, the matrix 
& in equation Eq. (|45|l may not be positive definite. 
This problem can be remedied somewhat by testing for 
positiveness using a Cholesky factorization, which takes 
negligible time compared with the para-unitary diagonal- 
ization. A better solution is to make a good initial guess 
to the values of {A„}. One can do this by first running 
the minimization with an almost uniform momentum dis- 
tribution (in case this is included as observables) or with 
(ajao) set equal to N to t/N, and then carrying out the 
subsequent program runs using the previously obtained 
{A„} as initial values. In the subsequent program runs, 
one then gradually lets the momentum distributions ap- 
proach the real distribution. A similar approach may 
be taken with the total number of particles, as making 
good guesses of the {A^} may be easier at low particle 
numbers. 

A trick we have made use of, which made guessing the 
{A^} easier, is scaling down of the measured data. Since 
we are not interested in the true many-body density op- 
erator for the system, but only the second moments of 
the ladder operators, we may as well make the calcu- 
lations for a smaller number of particles, provided the 
calculations are still precise, and then rescale the data 
afterwards. As an example, we shall show how to per- 
form this down scaling in the b-basis at t = 0. 

Let (fob') be the matrix of (unknown) second moments 
we wish to reconstruct. We will instead perform calcula- 
tions on the matrix (bb^)us, also having the character- 
istics of a matrix of second moments: 

(« , >«-«((« , >-( J J'S)) + ( / i r S) 

= «<»») + (1 - k) ( 7 J °) (Al) 
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where J/v is the NxN identity matrix and k is a real num- 
ber. When reconstructing the matrix (bb^)j)S we must 
modify the measurements, retaining gN to t and thereby 
the transformation matrices stf , 33 and . So, for in- 
stance, the down scaled spatial density distributions that 
should be used can be found from: 



DS 



^~ x 3§~ 



In 




-i\t 



Comparing the lower half part of the diagonals, we see 
that we should multiply the coordinate space density 
measurements by n and add to them ^2 v 2 . In 



a similar manner, measurements of other observables 
should be modified in the scaled down calculations. After 
the A„ has been found, the relation Eq. (|A1|) can easily 
be inverted, using the same point of time and basis as 
used in the original scaling. 



The two arguments n and q can both assume N values, 
in our case {-M, —M + 1, . . . , M}, like in Eq. The 
mapping is bijective so the Wigner function contains all 
the information in the one-body density matrix. Further- 
more, the wave number- and spatial densities are easily 
recoverable: 



M 



q=-M 
M 

(a\k q )a(k q )) = Y, W(n,q). 



n=—M 



Appendix B: WIGNER FUNCTIONS 

It is common practise to map the NxN one-body 
density matrix to the real Wigner function W(n, q) given 
on an N x N grid [l^ . 

1 M 

W(n,q) = — 22 p(f(n-y),f(n + y)) ■ 

y=-M 

exp (Amyn/N) (Bl) 
/(n - y) = mod (n - y, N) - M - 1. 



In the case of only one particle, the evolution of the quan- 
tum system is completely described by the Hamiltonian 
and the one-body density matrix. In our case, having 
many particles, we instead have to specify all second 
moments to know the evolution of these, see Eq. Q52p. 
Therefore, in this paper, the Wigner function is merely 
meant as an illustration of the one-body density matrix 
at a particular point of time. 
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The statement of the MAXENT procedure in this paper 
is very similar to the one given in f.x. [f| 
Going directly from Eq. Q, we have subtracted a con- 
stant, having no physical significance 
To follow [T3 we use the, in matrix algebra unusual, pre- 
fix "para-" instead of f.x. "pseudo-". 
Had we been treating fermion- instead of boson- 
operators, the matrix 2F would have been unitary. 
An bilinear, hermitian operator can always be brought to 
this form by adding a scalar constant which is determined 
by the diagonal of the matrix. 



